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Critical transitions in multi-stable systems have been discussed as models for a variety of phe¬ 
nomena ranging from the extinctions of species to socio-economic changes and climate transitions 
between ice-ages and warm-ages. From bifurcation theory we can expect certain critical transitions 
to be preceded by a decreased recovery from external perturbations. The consequences of this crit¬ 
ical slowing down have been observed as an increase in variance and autocorrelation prior to the 
transition. However especially in the presence of noise it is not clear, whether these changes in 
observation variables are statistically relevant such that they could be used as indicators for critical 
transitions. In this contribution we investigate the predictability of critical transitions in conceptual 
models. We study the quadratic integrate-and-fire model and the van der Pol model, under the 
influence of external noise. We focus especially on the statistical analysis of the success of predic¬ 
tions and the overall predictability of the system. The performance of different indicator variables 
turns out to be dependent on the specific model under study and the conditions of accessing it. 
Furthermore, we study the influence of the magnitude of transitions on the predictive performance. 


I. INTRODUCTION 

Evidences of abrupt drastic shifts, so-called critical 
transitions or tipping points , have been reported in very 
different complex systems including, but not limited to, 
ecosystems mm, medical conditions mm, financial mar¬ 
kets [5] and climate mm- When a critical transition is 
reached, the system undergoes a sudden, relatively rapid 
and often irreversible change. The idea that a wide va¬ 
riety of critical transitions (CTs) could be preceded by 
early-warning signs mi is summarized in the short re¬ 
view [10J . Early-warning signs are theoretically justified 
through understanding a system that is likely to undergo 
a critical transition as a dynamical system close to a bi¬ 
furcation point. In this case, perturbations from stable 
states decay more slowly in comparison to the situation 
far away from the bifurcation. This effect is well known as 
critical slowing down (CSD) 11- The consequences of a 
CSD can be monitored in several variables, that can serve 
to predict CTs. Examples for these predictors in stochas¬ 
tic systems are an increase of the variance PH and an 
increase in autocorrelation mm- This approach has 
been applied to models fTKIfTT] but also to experimen¬ 
tal data pans] and field data [20] . However, the above 
mentioned examples of critical transitions that have been 
predicted and observed only once mm do not allow us 
to access the quality of the predictors in any statistically 
relevant way. It is in fact common for indicators of ex¬ 
treme events EEH2II to be tested with standard measures 
for the quality of classifiers, such as skill scores [25, ;2S], 
contingency tables [30l or receiver operator characteris¬ 
tic curves (ROC curves) [30] ■ Additionally it is common 
to test for the dependence of the prediction success on 
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parameters related to the estimation of predictors, the 
prediction procedure and the events under study. Such 
parameters could be, for example, the length of the data 
record used to estimate the predictor, the lead time (time 
between issuing the forecast and the observation of the 
event) or the magnitude and relative frequency of the 
event under study [23]. Similar tests for indicators as¬ 
sociated with CSD are—apart from one study m — still 
missing. 


Although theoretical aspects of a relation between 
CSD and CTs can be formalized mathematically us¬ 
ing the theory of fast-slow stochastic dynamical systems 
Hi m, relatively little is known about the practical rel¬ 
evance of the indicators associated with CSD. There are 
several mathematical hypotheses in the theoretical re¬ 
sults, which may not be satisfied in practical applica¬ 
tions. In fact, real-life data sets are extremely difficult 
to analyze with respect to critical slowing down [ 1331435] . 
In the vicinity of many bifurcations, even relatively small 
perturbations can cause a drastic shift from a marginally 
stable state of the system to a contrasting state m, lead¬ 
ing to noise-induced transitions j36l 1'37| , which are more 
difficult to predict [38| as they are intermediate between 
a large deviation [3j5] and a bifurcation regime. 


In this contribution we evaluate the statistical rele¬ 
vance of predictor variables associated with CSD. In con¬ 
trast to single event predictions made for critical tran¬ 
sitions in experiments or large scale models, we study 
the success of predictors using data sets that contain 
10 6 CTs. To generate these data sets, we use two low- 
dimensional conceptual models, introduced in Sec. [TT] In 
Sec. |III[ we explain how we verify the occurrence of CTs 
in the generated time series. Details concerning the esti¬ 
mation of predictors associated with CSD are presented 
in Sec. IV The predictive power of these predictors is 
evaluated in Sec. [V] with a focus on varying the length of 
the observation used for estimating the predictors, and 
the lead time, i.e. the time between issuing the forecast 
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and the occurrence of the event. Additionally, we study 
in Secs. VI and VII whether the magnitude of the transi¬ 
tions has an effect of their predictability. We summarize 
our results in Sec. IVIIII 


II. CONCEPTUAL MODELS FOR CRITICAL 
TRANSITIONS 


One basic characterization of CTs is the separation 
of system into two distinct phases: (a) a slow drift of 
system variables without any large variations and (b) a 
fast transition phase with a very drastic sudden change of 
at least one of the system variables. Hence, it is natural 
to propose fast-slow dynamical systems uni as a basic 
building block to model these effects [32] ■ A deterministic 
planar fast-slow system of ordinary differential equations 
(ODEs) is given by 

= x = f(x,y), (1) 

^=y = g(x,y ), (2) 

where (x, y) £ R 2 , 0 < e <C 1 is a small parameter so that 
a; is a fast variable in comparison to the slow variable 
y. Although one may also study more general fast-slow 
systems with (x,y) £ R m x R", we shall restrict to the 
planar case as it already illustrates the effects we are 
interested in. The set 


C 0 := {(:r, y) £ M 2 : f(x, y) = 0} (3) 


is called the critical manifold; in fact, we are going to 
assume that Cq is a smooth manifold for the examples 
we consider. Considering the limit e —>- 0 in (jTJ) then 
leads to the slow subsystem on Cq given by 


dy 

dr 


V = 9 (h(y),y). 


(4) 


For e > 0, it is possible to rescale time in ®-§ via 
t := r/e and then consider e —► 0 to obtain the fast 
subsystem 


dx 

dt 


= x 


f(x,y), 


dy 

dt 


= y' = o, 


(5) 

( 6 ) 


which is a parametrized system with parameter y (or 
frozen slow variable). 

To incorporate the effect of noise, which may arise, e.g., 
via external forcing, finite-system size or internal small- 
scale fluctuations, one usually extends the system (|T|) and 
© by noise terms m and studies fast-slow stochastic 
ordinary differential equations (SODEs) given by 


X = -f(x,y) + — 7 = r ?i( r ), 
e V e 

(7) 

y = g(x,y) + <t 2 v2(t), 

( 8 ) 


where and 772 (t) are independent Gaussian white 

noises. 

There are different mechanisms that explain how large 
sudden jumps can be induced in fast-slow SODEs [55] , 
The two main classes are bifurcation-induced ( B-tipping ) 
or noise-induced ( N-tipping ) transitions. Regarding 
early-warning signs, it is well known that before B- 
tipping we expect a growth invariance and autocorrela¬ 
tion um- For example, the variance of the fast variable 
near the attracting part of the slow manifold approaching 
a fold point at (x c ,y c ) satisfies 

Var (x(t)) = O ^-^=L=^ , &sy\y c (9) 

as long as (x(t),y(t)) is outside a small e-dependent 
neighborhood near the fold point pb] i41] and suitable 
limits in e relative to the noise level are considered. A 
general variance increase up to a threshold and/or scal¬ 
ing laws such as ([9| can be used for predicting a jump 
at a fold point under suitable conditions on a sufficiently 
long data set and under the condition that the slow drift 
speed governed by e is s ufficiently larg e in comparison 
to the noise level cr := \/(oq) 2 + (cr 2 ) 2 ; for a fold, this 
condition is more precisely expressed by assuming that 
a = cr(e) depends upon e, 

0 < cr(e) as e —0 (10) 

However, if the condition fails and the opposite relation 
a \fe holds, we are in the N-tipping regime and the 
noise drives large jumps with high probability before the 
fold point governed by large deviation theory [39 . Purely 
noise-induced jumps are very difficult, if not impossible, 
to predict. 


A. The Quadratic Integrate-and-Fire (QIF) Model 

One natural model to start with is to restrict to a fold 
normal form with a reset after the jump at the fold; see 
Fig. dj In particular, if we let f(x 7 y) = y — x 2 , F(x, y) = 
1, g(x,y) = —1, and G(x,y) = 0 in ([7]) and Q, then we 
obtain a local fold normal form with noise acting on the 
slow variable 

x = -(y-x 2 ) + (11) 

e Ve 

y = -l. (12) 

We are also going to use a deterministic reset in a cer¬ 
tain phase-space region described below and refer to the 
resulting model as the quadratic integrate-and-ffre (QIF) 
model. It resembles the classical integrate-and-fire model 
1421 except that we use a quadratic function for the in¬ 
tegration phase before the firing phase. Near its bifur¬ 
cation point the QIF also resembles the theta model for 
excitable neurons [43l 144] . In the noiseless case, the value 
of y determines the number of equilibria of the fast flow 
of x. For a positive y 1 we have two equilibrium branches 
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FIG. 1. The dynamics of the QIF model illustrated by the 
phase portrait. The parameters (e,8,a i) = (0.02,0.5,0.2). 
The black curves represent the critical manifold C^ IF = 
{(x,y) : y = x 2 }, a solid line for the attracting part and 
dashed line for the repelling part. The fold point (0,0) is 
marked as a red (gray) circle. The gray dashed double ar¬ 
rows indicate the fast subsystem flow, i.e., when the system 
is not near C^ IF . A numerical solution is presented as a blue 
(gray) solid line, and its path of reset as a yellow (light gray) 
dot-dashed line. 


Cq = {x = ± s /y ) y > 0}, whose union together with 
the fold point at (0,0) is the critical manifold Cq. Note 
that Cq + is attracting, while Cq~ is repelling. When y is 
negative, the fast subsystem has no equilibria at all. In 
particular, a saddle-node (or fold) bifurcation occurs at 
y = 0. Therefore, the critical manifold of the QIF model 
is attracting in quadrant I (x > 0, y > 0) and repelling 
in quadrant II ( x < 0, y > 0) as is illustrated in Fig. [I] 

When 0 < e < 1 and starting from the point (1,1) 
and uniformly decreasing y , the trajectory of the solu¬ 
tion travels near the attracting critical manifold C^ + to¬ 
wards the fold point (0,0). Shortly before reaching the 
fold point, depending upon the noise level, the system 
may perform a noise-induced jump across the flatten¬ 
ing potential barrier between the stable and the unstable 
equilibria and arrive in quadrant II. In quadrant II, the 
repelling critical manifold drives the system further and 
further towards negative infinity in x. In our model, how¬ 
ever, the system is considered to be in an excited state, 
when the fast variable x is below a threshold —8 and reset 
to the initial condition (1,1). 

Numerical simulations of equations (111 and (12) us¬ 
ing the Euler-Maruyama method [45] generate time se¬ 
ries of N observations {£„} and {y n } at discrete time 
steps t n = to + nAt. Here, to denotes the initial time, 
n = 0,1,..., N — 1 is the index of each time step, and A 
is a constant time interval of numerical integration. 


As illustrated in Fig. [2j CTs can be observed in the 
time series of the fast variable x while the slow variable y 
acts as the slowly changing bifurcation parameter. Since 
the system is reset to the initial state (1,1) after x exceeds 



FIG. 2. Time series of the fast variable x and the slow variable 
y of the QIF model, simulated using the Euler-Maruyama 
method with parameter values (e, 8, ay) = (0.02,0.5,0.2). 


a certain threshold —<5, we can use the QIF model to 
generate an arbitrary amount of CTs, which we are going 
to investigate below from a statistical perspective. 

The stochastic QIF model may not only have direct 
relevance for many transitions in neuroscience [421 146j as 
a local normal form to model the subthreshold dynamics 
before spiking or bursting but the QIF model could also 
be viewed as useful for any applications with local fold 
dynamics and global resets. 


B. The van der Pol Model 

In addition to the purely local QIF model with resets, 
it is also natural to compare it to a model, where the 
resets are via a smooth global nonlinearity. The classi¬ 
cal example to study are van der Pol m (or FitzHugh- 
Nagumo [55|) relaxation oscillators [IS]- In particular, 
we consider f(x,y) = y - jgjX 2 {x + 8 ), F(x,y) = 1, 
g{x, y) = — | — x, G(x, y) = 0 and obtain a version of the 
van der Pol (vdP) system 

% = e {y~ w x2 ^ x+s ^) + ( 13 ) 

8 , , 

y = - ( 14 ) 

The precise choice of the form of the model will be mo¬ 
tivated in more detail below, particularly with respect 
to the parameter 8. When the external stimulus exceeds 
a certain threshold, the behavior of the system changes 
qualitatively from a stable fixed point to a limit cycle 
undergoing a Hopf bifurcation. 

The deterministic version of the model, i.e., <j\ = 0, has 
one fixed point (£fp,2/fp) = (—f> §f)> which is unstable 
under the assumptions 8 £ (—v/3, 0) and 0 < e C 1. 
A trajectory of the stochastic vdP (cti ^ 0) forms a 
noisy relaxation-oscillation-type periodic orbit involving 
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FIG. 3. The dynamics of the vdP model. The parameters 
(e, <5, <ti) = (0.02,0.5,0.1). Analogous to Fig. [l] the critical 
manifold (black lines, solid for the attracting part and dashed 
for the repelling part), the fold points at (—1<5,1) and (0,0) 
(red circles) and the numerical solution trajectory [blue (gray) 
solid line] are plotted in state space. The dashed double ar¬ 
rows indicate the orientation of the relaxations in the noisy 
case and the noise-induced transitions for the fast varia 



FIG. 4. Time series of the fast variable x and the slow variable 
y of the vdP model, simulated using the Euler-Maruyama 
method with parameter values (e, 5, ay) = (0.02, 0.5, 0.1). 


two rapid transitions and two slow drifts as is illustrated 
in Fig. [3j Since the critical manifold of the vdP model 
has two fold points at (— |<5, 1) and (0,0), the manifold 
is naturally split into three parts (left, middle and right) 


of a CT as the distance S for the jump at, respectively 
near, the lower fold. Analogous analysis applies for the 
second fold point as we have a fast subsystem at y = 1 
given by 


WdP = CvdP n <{ {x, y) : x < --S 


CZ P = C vd p D |(*,y) : --6 < x < 0 \ , 
CvdP = CvdP n {(x, y) : x > 0} . 


(15) 

(16) 
(17) 


By investigating the stability of the equilibria of of the 
fast variable x for a fixed y (in the e —> 0 limit), we see 
that C^ dP , C; dP are normally hyperbolic attracting parts 
of the critical manifold and CZ p is a normally hyperbolic 
repelling part of the critical manifold. Therefore, in the 
vicinity of C^ dP and CZp, the system is resilient to per¬ 
turbations for long times, i.e. it is able to relax quickly 
back to the critical manifold Furthermore, the sys¬ 
tem travels slowly along CZp and Cy dP , driven by the 
slow variable y. 

In quadrant I of the x-y plane, y is negative so the 
system moves downwards along CZp- Under conditions 
of the form (10) the system recovers from small random 
disturbances with high probability until it is close to the 
fold point (0,0). When near (0,0), the potential bar¬ 
rier between the stable equilibrium on C^ dP and on C^ dP 
flattens, so that it becomes more likely that the system 
jumps away from CZp to a region to the left of CZp and 
then it gets repelled to a region near the other attracting 
part Cy dP . Note very carefully that if we consider the 
fast subsystem precisely at y = 0 and for ay = 0, then 
we have 


x ' = {~^^ X + S ^) (18) 

so that the second fast subsystem equilibrium besides x = 
0 is given by x = — 6. So we may define the magnitude 


/ I i 27 2 

1 =|l -4ff X 


(x + S ) 


(19) 


with fast subsystem steady states at x = —|<5 and 
x = |(5, so the distance of the jump at the second fold 
is again of size S. This allows us to change <5 and as¬ 
sociate with it the size of both CTs during a relaxation 
oscillation. The only difference between the branches is 
that the speed of upward drift in the left branch is slower 
than the downward drift in the right branch, leading to 
a longer residence time for the left branch. 

Simulating Eqs. (131 and (fl4|) can again be used to 


generates time series of N observations {x n } and {y n } at 
discrete time steps t n ; see Fig. [4} In contrast to the QIF 
model no (potentially artificial) resets are necessary for 
the vdP model. 


III. DETECTING CRITICAL TRANSITIONS AS 
EXTREME EVENTS IN TIME SERIES 

The two types of CTs observed in the QIF and the 
vdP model, correspond to two types of extreme events: 
threshold-crossings and increments. Critical transitions 
in the QIF model can be categorized as threshold- 
crossings with the threshold — <5 (Fig. [5]). Relevant exam¬ 
ples for extreme events that consist in threshold crossings 
are river levels outreaching a flood warning threshold af¬ 
ter a heavy rain, or temperatures dropping below the 
freezing point. Increments, however are defined as the 
sudden change of the observed variable larger than a cer¬ 
tain value. Prominent examples for increments are surges 
in electric loads and drastic decreases of stock market in¬ 
dices. In the time series of the fast variable generated by 
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FIG. 5. Time series {x n } and tracer time series {x-n} of 
the QIF model with (e, <5, cri) = (0.02,0.5,0.2). In the up¬ 
per panel, the threshold for reset x = —5 is shown as a red 
solid line. Both time series are aligned such that the time 
points t n when the CTs occur (dashed lines) and the ones t m 
when the predictions are made (orange (light gray) impulses 
of the tracer time series) can be compared. The forecast lead 
time between f n and t m is set to k = 500 for demonstration. 
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FIG. 6. Time series {x„} and tracer time series {\>i} of 
the vdP model with (e, 8, ai) = (0.02,0.5,0.1). In the up¬ 
per panel, the transition phases are shaded with yellow (very 
light gray). Both time series are aligned such that the time 
points t n when the CTs occur (dashed lines) and the ones t n 
when the predictions are made (orange (light gray) impulses 
of the tracer time series) can be compared. The forecast lead 
time k is set to k = 1000 for demonstration. 


the vdP model CTs can be observed as increments with 
magnitude 5 (see Fig. [6). 

In order to monitor the occurrence of threshold cross¬ 
ings and increments we introduce a secondary time se¬ 
ries called the tracer time series {xn} with the time in¬ 
dex n = 0,1,... ,N and the same sampling interval A. 
The binary tracer variable \n stores information about 
whether an extreme event, i.e., here a CT happens at a 
future time point t n + kA. It is not used as an indicator 
variable for predicting future CTs, but only for marking 
the start time of CTs, as will be explained in detail in the 
following. The time kA, i.e., the time between a forecast 
made at t n and occurrence of an event at t n + kA, is 
sometimes called lead time of a prediction. 

Identifying threshold crossings in the QIF model and 
keeping track of their occurrence in a tracer time series 
is straightforward: Once the value of x drops below the 
threshold — S at t n , the occurrence of a threshold crossing 
is identified and the tracer variable \n{d,k) at t n — kA 
is set to 1: 


/ r 7 \ I 1 ■ If ^n-t-fc A d 

Xn{0,k) = \ Q 


if x . 


n+k — 


> -5. 


( 20 ) 


1 : if x n +k+M %n > d 

>0i 


Xn+(p— l)wid (^7 ^0 < 


with 


and 


A x n + k+M 


A x. 


—n-\-k 


0 


—n+fe ^ n 

■ x p—i < ^0 


A x™ +k > 


^ u o 


(transition upwards), 


or x n+ k+M - x n < -6 (21) 
A x n f +k+M < d 0 

—n-\-k —n+/c 


A Xn 


A x" +k < <?! 


—n-f-K a 
I ■" • 7 1 > a l 


(transition downwards); 
otherwise. 

n+k+M+Wf 


—n+k+M 

X f 


7r n + k — 


= A E 

1 


w id 


Wf 

J i=n+k+M 
n+k+pwid 

E < 

i=n+k+(p—l)wid. 


Identifying the onsets of the CTs in the vdP model re¬ 
quires us to distinguish among large fluctuations, meta 
stable states, and true CTs. We therefore introduce two 
auxiliary thresholds: 9q, the upper limit of the lower 
state, and 9 1 , the lower limit of the upper state. In addi¬ 
tion to the occurrence of a large change \x n +k+M — x n \ > 
S with respect to the current value x n , we also require 
that the finite-time average after the change x^ +k+AI re¬ 
mains above 9\ or below 9 0 (see Fig. [?|. In more detail, 
the tracer variable detecting the occurrence of a CT is 
defined as follows: 


Here, we use non-overlapping windows of Wid time 
steps and the values of the parameters, M , , Wf,9o 

and 6*i are set such that they capture the CTs in each 
specific data set. More specifically, choosing M to be 
slightly larger than the typical duration of the fast tran¬ 
sition phase for a CT and 6 to be larger than the size of 
typical fluctuations around each state allows us to cap¬ 
ture true CTs while excluding large fluctuations. Addi¬ 
tionally, testing for the future behavior after large incre¬ 
ments as summarized by x^ +k+M allows us to distinguish 
true CTs from transitions to meta-stable states. 
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FIG. 7. The detection method for positive and negative incre¬ 
ments in the vdP model. The time series {xn} is shown as blue 
(gray) lines and the finite time average prior to CTs as black 
steps. The finite-time average is used to locate the onset of the 
transitions by comparing it to the thresholds (dashed lines, do 
for the lower state and 6 1 for the upper state). The narrowed 
down transition phases (upwards-jumping in the left panel 
and downwards-jumping in the right panel) are colored yellow 
(light gray). The values of the parameters (M, Wid, uif,8 0 , 6 * 1 ) 
are chosen empirically as (300, 5, 20, —0.2, —0.1). 


The exact onset of a transition is determined as follows: 
We start from the point where a large change is detected, 
t n+ ky which should be somewhere before the transition, 
and calculate the finite-time average in non-overlapping 
windows. If in the p-th window starting from time t n +k 
the average x™ +k crosses the threshold, the onset of the 
transition is determined to be £„+*,+(p— (Eq. 211. 

The end of the transitions is determined analogously, i.e., 
starting from t n +k+M and calculating the finite time av¬ 
erage backwards until it crosses the other threshold. 


IV. PREDICTING CRITICAL TRANSITIONS 


Traditionally, one would think that the prediction of 
the future implies a precise modeling and understand¬ 
ing of the system under study, such as, e.g., in weather 
forecasting and climate modeling. However, whenever 
precise information about the system under study is un¬ 
known, or if the system is too complex to allow for a de¬ 
tailed numerically expensive modeling, then data-based 
predictions can be a suitable alternative US H 2 I- As 
more and more large data sets, machine learning meth¬ 
ods, and computation power become available [50] . data- 
based predictions can also be understand as a generic 
first choice that does not require human resources, but 
only the availability of sufficiently many data records of 
events. Apart from approaches that predict values of 
the future trajectories of systems edge], discrete events 
can be successfully predicted by observing relevant early- 


warning signals (precursors, predictors). Methods to 
identify relevant predictors and prediction margins con¬ 
sist in statistical inference ED [22H2H ESI ED ESI, artificial 
neural networks 15411551 or other pattern recognition and 
machine learning approaches M- Here we apply sta¬ 
tistical inference (naive Bayesian classifiers :53[ [57]) to 
predict CTs based on time series generated by the mod¬ 
els introduced in Sec. hd 


A. Estimating Indicators 


Knowing that CTs are preceded by a CSD, we can ex¬ 
pect all quantities that are affected by the slowing down 
to possess some predictive power. From a practical point 
of view, finite-time estimates of the standard deviation 
and the autocorrelation can be easily estimated from a 
time series. In more detail, we estimate the following 
indicators: 


PI a sliding window estimate of the variance ( v n ), 

P2 a sliding window estimate of the standard devia¬ 
tion (s„), 

P3 a sliding average of s denoted as (s) n and 


P4 a sliding window estimate of the lag-1 autocorrela¬ 
tion ( a n ) [551 . 

To mimic a realistic prediction scenario, we choose sliding 
windows of lengths wA such that the resulting indicator 
variable contains only information from the present data 
x n and from the past x n _ w+ i, x n _ w+ 2, x n _i, but not 
from the future. In other words, indicators at time t n , 
estimated from a sliding window between x n - w +\ and 
x n -i, are given by 


PI: 


P2: 

P3: 


P4: 


= v(x)t n = — Y] ( Xi-x tn ) 2 , 
w . _ 


i=n-w -\-1 


with Xt = ~ Y Xi 
n w . „ 

i=n-w -\-1 

s„ = s(x)t n = \Jv[x )t„; 

1 


S)n = 


S ( X )tn> 


i=n—uJavg + 1 


( 22 ) 


(23) 

(24) 


n —1 

F, (Xi - x tn ) (x i+1 - x tn ) 
a n ( 1) = i=ra ~ w+1 „ -■ (25) 

( X i~ X tr J 2 

i=n-w -\-1 


The indicators estimated from a sliding window estab¬ 
lish an indicator time series (Fig. [8]), which provides in¬ 
formation for prediction at each time step. The indica¬ 
tor values estimated from the time series including the 
critical transitions (shaded yellow in Fig. [8]) can be ex¬ 
ceptional, especially for the standard deviations. These 
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FIG. 8. Time series of tested indicator variables. From top 
to bottom, shown in the panels are the time series of x (vdP 
model with 5 = 0.5), of tracer variable x, an d of the four 
proposed indicators. The standard deviations and the auto¬ 
correlations are calculated in a sliding window with length 
w — 50 time steps. The window width for averaging standard 
deviations u> avg = 100 time steps. The tracer time series in 
the second panel is shown with a forecast lead time k = 500 
time steps. The time when predictions are made is marked 
with a dashed line in the four panels of indicators. The in¬ 
dicators within the intervals marked yellow (light gray) are 
excluded from the statistics for prediction. 


values have nothing to do with alarming the critical tran¬ 
sitions, therefore they are excluded from further analysis. 


B. Predictions 



Yn 

FIG. 9. An example for decision making via thresholds of 
CPDFs, as specified in Eq. ( |27| ). By gradually lowering the 
discrimination threshold and correspondingly extending the 
alarm volume, we can generate an ROC curve from a CPDF. 
The solid line with hollow squares represents a numerically de¬ 
termined CPDF. The discrimination threshold choosing high- 
performance indicators is shown as a dashed line and the cor¬ 
responding (disconnected) alarm volumes (filled rectangles) 
are indicated by black arrows. 


PDFs, which prevents over-fitting without discarding too 
much information. The most meaningful value of a in¬ 
dicator variable, the value most likely to be followed by 
an event, is the one that maximizes the CPDF. Relevant 
indicators should lead to non-flat CPDFs with one or sev¬ 
eral maxima. Apart from being sensitive to the occur¬ 
rence of the event, a meaningful indicator should also be 
specific, i.e., not occur by chance without being related 
to an event. One indication of specificity is that the max¬ 
imum of the CPDF does not coincide with a maximum 
of the marginal PDF. An alarm A n for CTs is raised at 
time t n if the value of the CPDF is above a discrimination 
threshold d , i.e., let V(d) = {tp n : P(Xn(i5) = 1 \ijj n ) > d} 
and consider 


A n bPn ,V(d)\ 


1, if ip n € V(d); 
0 , otherwise. 


(27) 


Gradually lowering the discrimination threshold d and 
determining the corresponding alarm volume V is illus¬ 
trated in Fig. [9] 


In order to distinguish between relevant and non- 
relevant values of indicator variables, we use a naive 
Bayesian classifier [53] i.e., a conditional probability dis¬ 
tribution function (CPDF), 


F (Xn(6) = 1| V>n), 


(26) 


where Xn{$) is the tracer variable as defined in Eqs. (20) 


and (21) and ip n denotes any indicator variable, such as, 
e.g., the ones defined in Eqs. (22)-(25). The CPDFs 


are obtained numerically from time series of indicator 
variables ip n by estimating the marginal P (if n ) an d the 
joint probability distribution function P(y n (5) = 1 ,ip n ) 
and then using Bayes’ theorem ||5S]. More precisely, the 
number of bins is fixed to 20 for both marginal and joint 


C. Quantifying Prediction Success 

A common method to evaluate the performance of a 
binary classifier is the receiver operating characteristic 
curve (ROC curve) [30 ED] • The idea of the ROC curve 
is to plot the rate of true positives r c (or hit rate ) versus 
the rate of false positives rf (or false alarm rate ) as the 
discrimination threshold of the classifier is varied. Each 
value of the threshold d corresponds to a single point in 
the ROC curve (r/,r c ) (see Fig. 10). The ideal classi¬ 
fier is expected to predict all events and issue no false 
alarms. Consequently, it should generate ROC coordi¬ 
nates close to the upper left corner, i.e., the point (0,1). 
The point (0,0) at the lower left corner of an ROC graph 
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FIG. 10. An example for an ROC curve (red dots connected 
by solid lines) and its summary indices. The AUC (area under 
the curve, shaded gray), and the KSD (Kolmogorov Smirnov 
distance, blue dotted line representing the distance between 
the upper left corner (0,1) and the ROC-curve) summarize 
the overall behavior of the ROC-curves. The two blue (gray) 
dashed lines connect the closest and the second closest points 
on the curve to the (0,1) point, which are needed in the calcu¬ 
lation of KSD. The corresponding likelihood ratio he., 

the slope of the ROC curve in the vicinity of (0, 0) is shown 
as a yellow dot-dashed line. The black dashed diagonal rep¬ 
resents an ROC-curve generated by random predictions. 


represents a value of d that is higher than any entry of 
the CPDF. A classifier whose predictive performance is 
better than random guesses corresponds to a point in 
the upper triangle with r c > 77 . The shape of the ROC 
curve characterizes the predictability of the positives and 
the quality of the applied classification strategy. One can 
show m that the slope of a proper |3DJ ROC curve must 
decrease as one moves up and to the right on the curve 
if decisions are made using a naive Bayesian classifier. 

To compare several ROC curves, scalar summary in¬ 
dices of ROC curves have been introduced. Common 
summary indices, such as the minimal distance to (0,1) 
(Kolmogorov Smirnov distance, KSD), the area under 
curve (AUC), and the slope of the ROC at (0,0) ( likeli¬ 
hood ratio) are illustrated in Fig. 10 In this contribution, 
the KSD is calculated as the distance between (0,1) and 
the line connecting the two closest points to (0,1) [see 


the dashed blue (gray) lines in Fig. 10 


V. COMPARING INDICATORS AND THE 
WAYS THEY ARE ESTIMATED 

Although theoretically postulated from the phe¬ 
nomenon of CSD m Isa, it is not apriori known how 
useful and specific indicators are connected to an in¬ 
crease of fluctuations in practical prediction scenarios. 
Additionally, one can also ask which indicator is most 
suitable for practical applications and how its predictive 
power depends on the way it is calculated and used and 
also on the events it aims to predict. Another particular 
aspect to test is the forecast horizon, i.e., how large the 
lead time k can be chosen without lowering the overall 
quality of predictions. Last but not least, we also test 
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FIG. 11. ROC curves generated by the indicators introduced 
in Sec. |IV A| used to predict CTs in the QIF model. In general, 
we observe an increase in prediction quality with decreasing 
lead time k and increasing window width w. 


how the quality of indicators depend on the width w of 
the time window that was used to estimate them from 
a given time series. In a realistic forecast scenario, i.e., 
when predicting only one CT from a relatively short data 
record, the width w might be limited by the total amount 
of available data. Working with models that create repet¬ 
itive CTs, and arbitrarily long time series, we can explore 
all useful variations of w. 


We therefore generate time series of 2 22 time steps from 
the QIF model and the vdP model, introduced in Sec. [Tl] 
and then estimate the respective indicator time series 
{ipn}- In more detail, we choose the additional parame¬ 
ter of the vdP model in Eq. (14) to be d = 0.5. The noise 
strength is chosen to be a\ = 0.2 for the QIF model 
and ( 7 1 = 0.1 for the vdP model. Choosing a smaller 
noise strength for the vdP model prevents the system 
from jumping back and forth frequently (flickering). For 
both models, the time separation parameter e is set to be 
0.02, which means that we are roughly in the intermedi¬ 
ate regime er = 0{^/e). To investigate the dependence 
of the predictability on the window width of the indi¬ 
cators, we estimate multiple indicator time series {7/7} 
for different window widths ranging from 10 to 1000. For 
the sliding average of s„, the window width for averaging 
u>stvg is set to be 100 time steps. In order to examine the 
influence of the lead time on the predictability, we gen¬ 
erate ROC curves for different k ranging from 1 to 1000 
time steps. When exploring the k-w space, the sum of the 
maximal values of k and w is restricted to the minimal 
interval between CTs. The ROC curves in Figs. m and 


12 illustrate the dependence of the prediction success on 
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FIG. 12. ROC curves generated by the indicators introduced 
in Sec. |IV A] used to predict CTs in the VdP model. In general 
we observe an increase in prediction quality with decreasing 
lead time k and increasing window width w. 


both variables, the forecast lead time k and the window 
width w used to estimate indicators. We observe that 
the prediction performance of all tested indicators and 
in both models decreases with increased forecast horizon, 
i.e., increased lead time. As observed for other predic¬ 
tion tasks, the smaller the difference between prediction 
and occurrence of the event, the greater the prediction 
success. The closer the system is to CTs, the stronger 
the early warning signals and the stronger is apparently 
the correlation of indicator and event. Indeed, a predic¬ 
tion extremely close to a CT (small k) may be very good; 
however, larger lead times would be highly desirable for 
applications. 

Additionally, we observe a dependence of the predic¬ 
tion quality on the way that indicator variables are es¬ 
timated. Using very small window widths leads to pre¬ 
dictions that are not much better than random. Increas¬ 
ing the window width w, and thus smoothing the in¬ 
dicator time series, all indicators show an increasingly 
strong ability to predict CTs. This increase is particu¬ 
larly dramatic for small lead times. In other words, close 
to the transitions, meaningful indicatory signals are more 
likely to be hidden within stochastic fluctuations. Conse¬ 
quently, increasing the window width can be an effective 
way to attenuate the influence of the fluctuations and to 
improve predictions. 


The dependence of the forecast performance on the 
leading time k and the window width w is clear in the 
heat maps of the ROC indices in the k — w parameter 
space (Figs. 13 and 14). For the QIF model (Fig. 131, the 
lag-1 autocorrelation a„(l) (left column) is the best indi- 
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FIG. 13. Heat maps of the ROC indices of the indicator 
o n (l) and s n in the k — w space for the QIF model. AUCs 
(Area under curves) are shown in the upper row and KSDs 
(Kolmogorov-Smirnov distances) in the lower row. 
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FIG. 14. Heat maps of the ROC indices of the indicator a n ( 1) 
and (s) n in the k — w space for the vdP model. AUCs are 
shown in the upper row and KSDs in the lower row. 


cator among the four proposed indicators: at the upper- 
left corner of the heat maps, where we make predictions 
close to the CTs but use a large window size, the AUC 
of a n ( 1) is the largest and the KSD is the smallest. The 
standard deviation s n performs slightly better than its 
sliding average (s) n and the variance v n . However, the 
performances of the three are qualitatively similar. 

For the vdP model, good predictions can be made even 
in a early stage, i.e. far from the CTs, if we use a large 
enough window to estimate the indicators (Fig. 14). In 
particular, the sliding average of the standard deviation 
(s) n (right column of Fig. 14) performs significantly bet- 
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FIG. 15. Magnitude dependence of the time interval between 
CTs. Dashed lines represent the average interval and solid 
lines the minimal interval between CTs. Green (gray) lines 
are for the vdP model and dark blue (black) ones for the QIF 
model. 


ter than the other indicators with smaller window widths 
for estimation (around 200 time steps): Even when the 
predictions are made far before the transitions they can 
be quite successful (AUC around 0.8 for k = 800). The 
performances of s n , v n are qualitatively similar to a n (l). 


VI. THE DEPENDENCE ON THE 
TRANSITION MAGNITUDE 


The magnitude of CTs generated by our models does 
depend on the model parameter S (see Fig. 15). In the 
QIF model, as introduced in Sec.|TT| this parameter coin¬ 
cides with the threshold on x for the reset of the system. 
In the vdP model, <5 determines the distance between the 
fold points and the corresponding attracting branches of 
the critical manifold to which the system jumps from a 
fold. In this section we address the question, whether 
larger transitions are proceeded by a more distinguished 
indicatory behavior and are therefore easier to predict. 
Previous studies on extreme events in stochastic pro¬ 
cesses have shown that prediction quality can sensitively 
depend on event magnitude and the decay of CPDFs in 
the limit of large events ;52]. We now test whether this 
idea can be linked to dynamical systems for CTs and 
their classical early-warning signs or not. As test data 
sets for each model, we use time series generated for 100 
different values of <5, i.e., starting from S = 0.5 to A = 1.5 
with increments of 0.01. To detect CTs in these time 
series (see Sec. IV AI the detection methods have to be 
adapted according to each value of 5. For the QIF model, 
the threshold for identifying CTs is simply set to be —5. 
For the vdP model, the parameters and 6 \ for locat¬ 
ing the onset and the end of the transitions were (after 
empirically testing) set to 


9 0 = 6* 0 - 0.6 • (5 - 0.5) (28) 

e 1 = 01 -0.1 ■ ((5-0.5), (29) 



0.2 0.4 0.6 0.8 1 

False positive rate 


k = 500 
. w = 10 

/(b) 

- 



6 = 0.5 —I— - 

*// 

5 = 0.7 - - 

6 = 0.9 •••■•*• 

Jr 

8 = 1.1 — • — 

8 = 1.3 • — ^ ■ • 

pi- 5 = 1.5 

jjr random guess —- 


0.2 0.4 0.6 0.8 

False positive rate 



0.2 0.4 0.6 0.8 1 

False positive rate 


8 = 0.5 ■ 
8=0.7 
8=0.9 
8 = 1.1 
8=1.3 -5 
8= 1.5 

random guess — - 
0.2 ~0A 06 08 

False positive rate 


FIG. 16. Magnitude dependence of the ROC curves (a)-(d) of 
Sn in the QIF model for different combinations of lead times k 
and window widths w. For small lead times we observe an in¬ 
crease in prediction quality with an increasing CT magnitude 
5. 


with 0 q = — 0.2 and 0* = —0.1 being the thresholds for 
S = 0.5 [see also Eq. (21) and Fig. [T]. 

We focus on the magnitude dependence of the standard 
deviation and evaluate the quality of predictions using 
ROC-curves for four representative combinations of small 
and large lead time and window widths. 

In total, we explore the three-dimensional parameter 
space of k, w, and <5, which is a key practical contribution 
and has been missing in previous work on CTs. The 
maximal values of k and w are set to be 500, so that 
the sum of k and w does not exceed the minimal interval 
between transitions. We also observed that the minimal 
interval between transitions decreases by about a half 
with increasing S for the vdP model (Fig.[l5|). This effect 
can be understood as a consequence of an increase in the 
velocity of the slow variable y at the left branch y £ 
[|<5, | J], which drives the system towards the fold point 
(—1<5,1) [see also Eqs. (13) and (|T4|)]. 


A. Results for the QIF Model 

Starting from time series of CTs and the indicator vari¬ 
able s n generated by the QIF model, we predict CTs and 
evaluate the quality of the predictions using ROC curves 
in Fig. [16] and the related indices (AUC, KSD, and like¬ 
lihood ratio) in Fig. [l7] If the prediction is made close 
to the transition (small lead times), we observe a clear 
positive magnitude dependence, i.e., larger events can be 




















11 



5 8 


FIG. 17. Magnitude dependence of the predictability indices 
of the QIF model. ROC-indices, AUC and KSD, are plotted in 
(a) and (b), and the best indicator and the likelihood ratio 
evaluated at the best indicator in (c) and (d). The indices for 
k = 1, w = 10 are shown with solid lines (—), for k = 1, w = 

500 with dot-dashed lines (-), for k = 500, w =10 with 

dashed lines (-) and for k = 500, w = 500 with dotted 

lines (■ • •). For small lead times we observe an increase in 
prediction quality with an increasing CT magnitude 5. 


better predicted as shown in Figs. |T0[a) and |T0[c). For 
larger lead times the prediction quality decreases and the 
dependence on the magnitude vanishes completely. The 
ROC indices in Fig. |17| nicely summarize these results by 
showing an increase in AUC and likelihood ratio and a 
decrease in KSD with increasing 8 for small lead times. 
For larger lead times the indices do not show any clear 
trend and they are, up to small fluctuations, constant. 

One can understand the improvement of the predic¬ 
tion quality close to the transitions (i.e., with small k) as 
a consequence of the accelerating drifting of x towards 
negative infinity (see Fig. [I]) before the CT happens. As 
can be seen from the model’s dynamics [see Eq. 0 and 
Fig.g the drifting speed of the fast variable x rises as the 
system has passed the fold point (0,0). The larger the 
reset-threshold 5 the longer this faster drifting continues, 
before the system is reset to (1,1). As a consequence, 
the sliding s n at k time steps before the reset (i.e., the 
occurrence of a CT) becomes larger and larger, due to 
the accelerating drifting of x. In fact, the values of the 
maximum CPDF indicators s* increase with increasing 
threshold 8 as can be seen in Fig. 17 )c). 

Since these large indicatory s n are easier to distinguish 
from random fluctuations in the s n , also the likelihood 
ratio estimated for s n s, which are close to the transition 
(small lead times), increase with 5 and are up to three 
orders of magnitude larger close to the event [as shown 
in Fig. [Tt]T 1)] than likelihood ratios for large lead times. 
In contrast to this, s n s further away from the CT are not 
affected by the increase of the reset-threshold 8 but follow 
the (^-independent) dynamics, given by Eqs. (11) and 
(12). Consequently, we observe a sensitive dependence 




FIG. 18. Magnitude dependence of the ROC curves (a)-(d) 
of s„ in the vdP model, for different combinations of lead 
times k and window widths w. For large lead times and large 
window widths (d) and small lead times and small window 
widths (a) we observe a mild decrease in prediction quality 
with an increasing CT magnitude 8. 


on the event magnitude only for small lead times. 


B. Results for the vdP Model 


Testing for the dependence of the prediction quality on 
the event magnitude in the vdP model, we find qualita¬ 
tively different results. Instead of improving with increas¬ 
ing transition magnitude 8 , the quality of predictions re¬ 
mains the same or slightly decreases (negative magnitude 
dependence) as illustrated in Fig. 18 The dependence on 
8 is especially pronounced for the combinations of small 
k and small w and large k and large w. For the combina¬ 
tions of small k and large w or large k and small w, there 
is no clearly visible dependence on the transition mag¬ 
nitude. Two of the ROC indices in Fig. 19 (AUC and 
KSD) support these findings. However, the likelihood 
ratio shows a decrease for small window widths and no 
clear trend for large window widths. This deviation can 
be understood as a consequence of the definition of the 
likelihood ratio. Since it is evaluated at the value of the 
maximum CPDF indicator s* it reflects the slope of the 
ROC curve in the region of small false positive rates, i.e., 
close to the origin (0,0). In total we can conclude that 
the predictability of the vdP model is either independent 
or decreases with increasing magnitude of the transition. 
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FIG. 19. Magnitude dependence of the predictability indices 
of the vdP model. ROC-indices, AUC and KSD, are shown in 
(a) and (b), and the best indicators s„ and the likelihood ratio 
evaluated at the best indicator in (c) and (d). The indices for 
k = 1, w = 10 are shown with solid lines (—), for k = 1, w = 

500 with dot-dashed lines ( - ), for k = 500, w =10 with 

dashed lines ( - ), and for k = 500, w = 500 with dotted 

lines (•■■). For large lead times and large window widths 
and small lead times and small window widths we observe 
a mild decrease in prediction quality with an increasing CT 
magnitude 5 in (a) and (b). 


VII. TOWARDS A MORE DETAILED 
UNDERSTANDING OF THE MAGNITUDE 
DEPENDENCE 






FIG. 20. Magnitude dependence of the QIF model according 


to the condition c(ip n ,S) [see Eq. (301] for four representative 


combinations of k and w. The values of c(p n ,S) are repre¬ 
sented as grayscale-coded dots with a black border as long as 
they can be evaluated. The location of the optimal indicator 
are marked as a “x” if the value of c(ipp(8 ), 5) is available, 
otherwise and it is marked as a black . Note that c(ipn,5i) 
can only be evaluated, if PDFs, CPDFs and their derivatives 
could be estimated with sufficient accuracy for all values of 
ip n , Si, and 5i+i, where the index i enumerates discrete steps 
in 5. 


While the likelihood ratio Ais a useful summary 
index for analytically determined ROC curves, it is less 
reliable for numerically determined ROC curves, which 
may not be smooth. Nevertheless, A[^/>*,Xn(^)] can still 
be used to investigate the magnitude dependences of the 
best indicators. Reformulating the partial derivative of 
the likelihood ratio with respect to <5, one obtains the 
condition [55] 


cWn,6) := d 5 In [L^ n ,S)\ - 1 - d s ln[P(<T)]. 

(30) 

Here, L(ip n ,8) = P(Xn+fc(£) = and P(S) denotes 

the total probability of finding transitions which are 
larger or equal than S. A positive c at ip n = ip^ indicates 
that the best precursors’ prediction quality improves with 
larger event magnitudes, and a negative one indicates 
the opposite. Numerical estimates of c = c(ip nj 5) are 
presented for discrete values of %p n and S in the ifi n — 6 
plane (Figs. 20 and 21). The derivatives ds ln[L(ip n , 5)] 
and ds ln[P(o]j are numerically estimated using PDFs 
and CPDFs for neighboring values of Si and d,;+i ■ Addi¬ 
tionally, all numerical estimates of distribution functions 
and their derivatives are smoothed using Savitzky-Golay 
filtering |64j of fourth order with 20 supporting points 
to the left and to the right. Note that the condition in 
Eq. (30) can be only evaluated, if all numerical estimates 


of distribution functions can be obtained with sufficient 
accuracy for all values of ip n and each neighboring pair of 
Si and Blank (green/medium gray) areas in Figs. 20 
and [2l] indicate that this is not always possible, given the 
number of events for specific values of %p n and S. Values 
of c(ip n ,5) are represented by grayscale-coded dots with 
a black border, locations of optimal indicators m 

ip*(5) := argmaxP(x„ = 1| ip n ,8), (31) 

4>n 

are indicated by crosses: If the value of c[ip* (5), <5] can 
be evaluated, then it is shown as a grayscale-coded “x”; 
if not, the location of ^/>*(<5) is shown as a black 
The results for the QIF model ( |20| are mostly consistent 
with the previous findings presented in ROC plots (see 
Fig. 16). If the lead time is small [k = 1, see Figs. 20 
(a) and 20 (c)], then c{ip n ,8) is positive for relevant indi¬ 
cators. For larger lead times (k = 500) there are almost 
equally many positive and negative values of c{ip ^, 5) [see 
Fig. [20] (b) and [20] (d)]. 

In contrast to this, the results for the vdP model 


(Fig. 21) indicate different magnitude dependences than 
the ROC curves in Fig. EL For small window widths 
[«> = 10, see Fig. 1211(a) and 21 (b)] no clear preference of 
the sign of c(ip*,o )'can be seen while the ROC curve for 
k = 1 ,w = 10 shows a negative magnitude dependence. 
If the window width is large [w = 500, see Fig. 21 (c) and 
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FIG. 21. Magnitude dependence of the vdP model according 
to the condition c(i/j n ,S) [see Eq. (301] for four representative 
combinations of k and w. Color coding as in Fig. [20] 


21 (d)], then it is not possible to evaluate Eq. ( |30| for op¬ 
timal indicators i/i*, indicated as black crosses. However, 
in the regions where it can be evaluated the c(tp n ,6) is 
negative. 

In general, it is a drawback of the condition in Eq. ( [30| 
that it is evaluated only for one specific indicator value, 
since in applications alarm volumes containing ranges of 
useful indicator values are relevant [see, e.g., Eq. (27)]. 
Additionally, we saw that the maximum CPDF indicator 
can also vary with the event magnitude 5. Therefore, we 
propose a quantity that summarizes Eq. (30) for a couple 


of suitable indicators. This combined alarm volume V* 
contains a range of indicator values comprising all best 
indicators ?/i*(8) [see Eq. ©] for all given range of 5, 


i.e., 


V* := [min(5), max ip^{5)\. 

O O 


(32) 


Analogously to Eq. (30), a condition involving the com- 
i be obtained using the dei 
)od ratio, i.e., 

'/y. dV>„P(V’n|Xn(<*) = 1)' 


bined alarm volume can be obtained using the derivative 
of the combined likelihood ratio, i.e., 


d s A* (V*, X n(6)) = d s 


Jy. dV>„P(V’n|Xn(<5) = 0) 


.(33) 


Since we are interested only in the sign of this expression, 
we can formulate the condition 


c v{V* : = Wi — hdglo, 


(34) 


with Ii = 
and Iq = 


'V* 


dlpnP (^nlXniS) = 1 ) , 
dljjnP (lpn\Xn(S) = 0) . 


IV 


We evaluate the integrals Ii(S), Io(S) numerically and 
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FIG. 22. Magnitude dependence of the QIF model using com¬ 
bined alarm volumes c* := cv(V*,x(S)) [Eq. (341] for four 
representative values of k and w. 
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FIG. 23. Magnitude dependence of the vdP model using com¬ 
bined alarm volumes c* := cv(V* ,x(d)) [Eq. ( |34| )] for four 
representative values of k and w. 
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estimate their derivatives by using the Savitzky-Golay fil¬ 
tering as described above. It is worth mentioning that the 
combined likelihood ratio A* (E*,x(<5)) is not the slope 
of the ROC curve, but an approximation of it, since the 
alarm volume for an ROC curve is defined through the 
discrimination threshold d [see Eq. (27)], not the com¬ 
bined alarm volume. Plotting cy as a function of S (see 
Figs. 22 and 23), we find that the sign of cy corresponds 
well to the magnitude dependence of the ROC curves in 
Figs. [T6] and [18] mostly positive for k = 1 in the QIF 
model and mostly negative for the vdP model. 


VIII. CONCLUSIONS 


CTs are abundant in nature and society and can have 
drastic consequences. In contrast to forecasting future 
values of systems, predicting CTs on the basis of early- 
warning signals can be understood as a classification 
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task. From theoretical considerations one can expect that 
many classes of CTs are preceded by a vulnerability to 
external perturbations, i.e., by CSD UM- As a conse¬ 
quence of CSD, several indicator variables such as lag-1 
autocorrelation and variance were proposed EM- In ap¬ 
plications, these indicator variables have been observed 
to show a specific behavior previous to the observation of 
one single CT Ml- However, the general quality of these 
indicators, their performance in the presence of noise and 
in the limit of short observation time series has so far not 
been accessed in a statistically significant way. In this 
contribution we have evaluated the quality of indicators 
derived from CT by using long time series of many CTs, 
generated by models of two different fast-slow systems. 
We have quantified the prediction success by computing 
ROC curves and related summary indices. 

We found mechanisms influencing the prediction suc¬ 
cess to be either sensitive to the specific model or model- 
independent. Effects that were observed in both mod¬ 
els are the dependence of the length of the data record 
and the forecast horizon: As common in forecasts, the 
prediction performance of all tested indicators and in 
both models decreases with increased forecast horizon. 
Increasing the window width, and thus smoothing the 
indicator time series, all indicators show an increasingly 
strong ability to predict CTs. Consequently, increasing 
the amount of available training data can be an effective 
way to improve predictions. Surprisingly, we find that 
the performance of indicator variables is dependent on 


the specific model under study and the conditions of ac¬ 
cessing it. For the QIF model, the lag-1 autocorrelation 
produces the best results in our numerical study in the 
limit of small lead times and large window sizes for esti¬ 
mation. However, for the vdP model, the sliding average 
of the standard deviation performs significantly better 
than other indicators in the limit of low data availability 
and long lead times, with a remarkable AUC score of 0.8 
or higher. 

Also the role of the transition magnitude depends sen¬ 
sitively on the model under study. For the QIF model, 
we observed a positive magnitude dependence for short 
forecast horizons, while larger transitions cannot be pre¬ 
dicted more easily if the forecast horizon is large. For CTs 
generated by the vdP model, the quality of predictions 
can even decrease with an increasing transition magni¬ 
tude. Consequently, high-impact events can not always 
be expected to be easier to predict, although results for 
short forecasts might indicate this. 
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